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Abstract 

The force required to pull 7t and 40t polyhedral boulders from the surface of an 
asteroid is simulated using the discrete element method considering the effects of 
microgravity, regolith cohesion and boulder acceleration. The connection between 
particle surface energy and regolith cohesion is estimated by simulating a cohesion 
sample tearing test. An optimal constant acceleration is found where the peak net 
force from inertia and cohesion is a minimum. Peak pulling forces can be further 
reduced by using linear and quadratic acceleration functions with up to a 40% 
reduction in force for quadratic acceleration. 

Introduction 

Based on the concept of redirecting a Near Earth Asteroid (NEA) to the Earth-Moon 
vicinity (Brophy, Culick, Friedman, Allen, Baughman, & Bellerose, 2012), NASA is 
developing the Asteroid Redirect Mission (ARM) to identify, capture and redirect an 
asteroid to a safe lunar orbit for future exploration. Two concepts are proposed. The 
first concept is to capture and redirect an entire small asteroid and bring it to lunar 
orbit for future study by astronauts and robotic missions. The second concept is to 
retrieve a boulder from a larger asteroid and return it to the same lunar orbit. 

The boulder retrieval operation contains many uncertainties. One of the problems that 
arises is the need to identify the forces required to separate the boulder from the 
surface of the asteroid. The magnitude of these extraction forces can be affected by 
the size and shape of the boulder, the strength of the regolith and the depth to which it 
covers the boulder, the acceleration imparted to the boulder during extraction, and the 
mass forces of gravitational attraction and apparent centrifugal force due to asteroid 
rotation. Regolith strength is affected by its packing density, grain size and shape 
distribution, grain contact friction, degree of grain interlocking, gravitation attraction 
of each particle with the mean field gravitational force of the rest of the asteroid, 
electrostatic forces, solar radiation pressure, and van der Waals cohesive forces 
between regolith grains and between the grains and boulder. The strength can also be 
affected by cold welding and the presence of water or ice in the asteroid as recent 
observations of water vapors around Ceres (Kuppers, 2014) may indicate. Van der 
Waals cohesive forces are considered to be a dominant contributor to regolith 
strength, competing with regolith weight, but greater than electrostatic and solar 
radiation pressure forces (Scheeres, Hartzell, Sanchez, & Swift, 2010). 


The goal of this work is to use the COUPi (Controlled Objects Unbound Particles 
interaction) discrete element method (DEM) model (Kulchitsky, Nye, & Johnson, 
2013; Nye, Kulchitsky, & Johnson, 2014) to perform simulations of boulder retrieval 
in the asteroid surface reference frame and calculate the reaction forces. 

Contact Model Description 

In the COUPi DEM model, the contact forces are represented by normal elastic forces, 
adhesive forces, tangential forces and normal and tangential damping forces. The 
forces are computed using a modified Hertz-Mindlin model. The detailed explanation 
of COUPi DEM contact physics is described in (Nye, Kulchitsky, & Johnson, 2014). 
The normal force component depends on the overlap between particles (6). The 
elastic component of the contact force can be found as follows: 
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where the effective Young’s modulus and effective interaction radius r* are 
defined as follows: 

r, = (1/ri + l/r^y^ , E, = ((1 - vl)/E^ + (1 - vD/E^Y^ 

where and r 2 are the radii of contacting particles curvature at the contact point, 
and V 2 are Possion’s ratios, and and E 2 are Young’s moduli of contacting 
materials. The normal force damping component is used from the extended Hertz 
theory (Kuwabara & Kono, 1987): 

Eni = kfiiS^^^dSIdt, 

where the constant damping coefficient /c„j nonlinearly depends on the coefficient of 
restitution (Muller & Poschel, 2011). 

The adhesive normal component of the contact force is calculated using Johnson, 
Kendall, and Roberts (JKR) model (Johnson, Kendall, & Roberts, 1971; Johnson, 
2003) with linearization proposed by (Brilliantov & Poschel, 2004; Brilliantov, 
Albers, Spahn, & Poschel, 2007): 

= (4V^ /9) 7l/2£-U2^^3/4^3/4^ 
where y is surface energy value for interacting materials. 

The tangential force value at each time step t + At is calculated incrementally 
from the force at the previous time step t similar to (Zhao, Nezami, Hashash, & 
Ghaboussi, 2006; Nye, Kulchitsky, & Johnson, 2014) using the following expression: 

0t+At = ft - (Rt ■ nOn^ - kyy^^At, 

where jU is the Coloumb coefficient of friction between two particles, k-^ is the 
tangential stiffness, is the normal vector to the contact plane at time t, and is 
the relative velocity of the particle at the point of contact. 




Figure 1. Tri-spherical particles used in simulations. 


Particle Shape 

To simulate the regolith, we used a monodisperse system of partieles ereated with 3 
equal spheres symmetrieally conneeted together sueh that they have only one 
eommon point as shown in Figure 1. 

Cohesion Scaling 

Introduction 

The eohesive foree at the eontact in the JKR model is defined by the surface energy 
parameter y and the effective curvature radius r* of the particles at the point of contact 
and depends on the type of the contact and other material and contact properties. It 
can be shown (Johnson, 2003) that the maximum attracting force at the contact in 
JKR model is proportional to both y and r* and does not depend on other material 
properties: 

Any pulling force stronger than this will result in the failure of a static contact. Thus, 
Fjjjax is the strength of the cohesive force between the particles. Strength Fj^^x is a 
local parameter that determines the strength of a single contact. This parameter needs 
to be linked to the macroscopic regolith cohesive strength defined as the maximum 
pulling stress a regolith sample can have at any point (Li, Zeng, & Wilkinson, 2013). 
The stress in granular material is defined on a surface within the granular material 
that intersects a large enough number of particles. 

DEM simulations are performed with particles that are significantly larger than actual 
regolith particles due to computational limitations. As particles become larger the 
cohesive forces become smaller according to Eq. ( 1 ). Thus the surface energy y 
needs to be scaled up to accommodate this difference in size. For small values of 
surface energy, the standard triaxial compression test might not be sensitive enough 
(Modenese, Utili, & Houlsby, 2012). To link the cohesion stress with y andr,, we 
designed a cohesion gravity test as described below. 
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Figure 2. Cohesion gravity test. Left panel: initial setup, Right panel: final stage. 

Cohesion test description 

The main goal for this test is to measure the separation cohesion stress in granular 
materials depending on the surface energy and particle size. The test uses an 
increasing gravitational field to ramp up a uniaxial internal tensile stress in the 
regolith to induce failure. We keep the packing density the same for all the tests. The 
cohesion test consists of the following stages: 

1. Filling the particles. The goal for this stage is to fill a box with 27 groups of 
particles such that each group forms a connected rectangular specimen as 
shown in Figure 2 (left panel, each group shown with its own color, the 
specimen is shown upside down). To do this, 3 layers of identical tri-spherical 
particles with 9 groups in each layer are filled one by one. To fill a layer, the 
particles are put into the grid above the previous layer or the bottom of the 
box so they do not overlap. Then, they fall under the gravity 10 m/s^. At this 
stage there is no cohesion set between the particles, e.g. y = 0. The Coulomb 
friction is chosen to control the packing ratio for the particles. 

2. Cohesion application. The goal of this stage is to equilibrate the particles 
with surface energy set to the desired values. The surface energy is a property 
of the materials interacting and can be set up for each pair of materials in the 
system. We use the following rules for setting the surface energy: 

a. The two upper layers consist of particles with the desired surface 
energy y. There is neither friction nor adhesion set between those 
layers and the box walls to avoid their interaction. 

b. The 9 groups of the lowest layer (shown on top in Figure 2) have 10 
times higher surface energy between its particles. It has normal surface 
energy y between its particles and the particles from other layers. It has 
10 times higher cohesion between the layer and the wall of the box. 
The meaning of this set up is explained below. 

The material equilibrates for a few seconds to adjust to the new surface 
energy values. 

3. Gravity relaxation. The goal for this stage is to remove the gravitational 
forces. Due to the potential energy stored in each contact, the gravity cannot 
be removed immediately. Instead, a linear gravity relaxation from initial value 
oflOm/s^toO within 10 s is applied. After this the specimen is left for 20 s 
to equilibrate in 0 gravity conditions. 



Figure 3. Cohesion stress as a function of surface energy for different radii. 

4. Gravity pulling. At this stage a negative gravity that is linearly increasing 
with time is applied until the sample breaks as shown in Figure 2 (right panel). 
As particles on the lowest level have significantly stronger bonds, the 
breakage occurs between layer 1 and layer 2 where the maximum hydrostatic 
pulling stress occurs. Moreover, the particles that are strongly connected to 
the box do not separate from the box. The net force vector between central 
groups of particles between layer 1 and layer 2 is recorded as a function of 
time. The pulling stress between those groups is calculated by dividing the net 
force acting from the particles of one group acting on the particles of the other 
group by the cross-sectional surface area separating those groups. The 
maximum pulling stress is recorded as a cohesion stress for this specimen. 
The groups are chosen in the center to avoid any possible boundary effects as 
shown in Figure 2. 

Cohesion test results 

We tested the particles from 0.025 m to 0.05 m to find how cohesion stress depends 
on the surface energy and the radius. The results are shown in Figure 3. We found 
that for the range of cohesion stress we are interested in and for a packing density 
between 0.57 and 0.58, a good approximation for cohesion stress can be given by the 
following formula: 

P = Vo(y/Yof'^ 

where Po, 7o» ^^e some constants as shown in Figure 3. 

Simulation Setup 

The simulation consists of the following stages. 

1. Filling the particles. The goal for this stage is to fill a box with particles. The 
particles are put into the grid above the bottom of the box so they do not 
overlap. They are then allowed to fall due to gravity at 10 m/s^. At this stage 
there is no cohesion set between the particles, e.g. y = 0 . The Coulomb 
friction is chosen to control the packing ratio for the particles as in the 


Table 1. Material and contact properties and simulation parameters used in the tests. 
The material stiffness was reduced to increase the time step for the computations. 


Material/Contact/Simulation Property 

Value 

Material density, p, kg/m^ 

2875.0 

Packing density 

0.57 (approximately) 

Stiffness, G, GPa 

0.1 

Poisson coefficient, v 

0.25 

Interparticle friction ratio, p 

0.3 

Coefficient of restitution (at 1 m/s), Cr 

0.1 

Tangential contact stiffness, kj 

0.7 

Boulder mass 

7t and 40t 

Boulder equivalent sphere diameter 

1.67m and 2.98m 

Number of particles for cohesion tests 

38K, 81K,and302K 

Number of particles for boulder 7t test 

27,380 (82K spheres) 

Number of particles for boulder 40t test 

176K (528K spheres) 


cohesion test. The polyhedral boulder of a realistic shape is planted in the 
upper layers and is deposited together with the particles. 

2. Cohesion application. The goal of this stage is to equilibrate the particles 
with surface energy set to the desired values. 

3. Gravity adjustment. The goal for this stage is to relax the gravity forces to 
the forces estimated on the asteroid. Similar to the cohesion test, the linear 
gravity relaxation from initial value oflOm/s^ to 2.5-10“^ m/s^ within 
10 s is applied. After this the specimen is left for 20 s to equilibrate in new 
gravity conditions. 

4. Pulling the boulder. The boulder is pulled with controlled vertical velocity. 
The forces acting on the boulder from the regolith are recorded. The vertical 
velocity can be an arbitrary function of time and were set to grow linearly, 
quadratically or cubically to a particular constant value. 

The material, contact and simulation parameters are shown in Table 1. 

Tests 

The tests were performed for two sample boulders with mass 7t with equivalent 

sphere diameter of 1.67 m and 40t with equivalent sphere diameter of 2.98 m based 



Figure 4. Pulling the boulder visualization. 


on estimations of the boulder sizes suitable for ARM mission (Mazanek & Merrill, 
2013). The goal was to establish the maximum pulling force depending on the 
estimated cohesion stress of the regolith as well as the time dependence of the pulling 
force for different methods of pulling. The boulder was planted shallow with less than 
30% being under the surface. Then, vertical acceleration was applied until the boulder 
reached the velocity of 0.2 m/s. The force vector acting on the boulder from the 
particles was recorded. The absolute value of the force vector is considered to be the 
force applied to the boulder due to cohesion in the regolith. The test configuration is 
shown in Figure 4. 

Table 2. Simulations of the 40t (40,000 kg) boulder pull. 


Test 

Y, J/m^ 

Pc Pa 

tfi s 

a, m/s^ 

Vf, m/s 

FcN 

F’t,N 

1 

165 

25 

4 

0.05 

0.2 

287.22 

2287.22 

2 

915 

350 

4 

0.05 

0.2 

1676.00 

3676.00 

Table 3. Simulations of the 7t (6974 kg) boulder pull. Notation: /7c is estimated 
cohesion stress, t/is total pulling time, a is acceleration, Vf is final velocity of the 
boulder, Fc is maximum force from the cohesion. Ft is the maximum total force 

Test 

Y, J/m^ 

Pc Pa 

tfi s 

a, m/s^ 

Vf, m/s 

FcN 

F’t,N 

1 

165 

25 

8 

0.025 

0.2 

48.48 

222.83 

2 

460 

125 

8 

0.025 

0.2 

237.45 

411.80 

3 

915 

350 

8 

0.025 

0.2 

958.75 

1,122.1 

4 

915 

350 

4 

0.05 

0.2 

615.62 

964.32 

5 

915 

350 

2 

0.1 

0.2 

691.22 

1,388.62 

6 

915 

350 

8 

0.25t 

0.2 

824.53 

928.12 

7 

915 

350 

8 

6.25 ■ 10“3tV2 

0.2 

524.33 

555.00 

8 

2600 

1650 

8 

0.025 

0.2 

2330.42 

2504.77 


The primary set of tests was performed with constant vertical acceleration. A Linear 
and quadratic vertical acceleration pull functions were also applied for 7t boulder. 
The results of the simulations are represented in Table 2 and Table 3. The cohesion 
strength was varied from 100 Pa to 2000 Pa to cover a wide range of potential 
extraction forces. 

Discussion 

The simulation runs for the 7t boulder with cohesion strength of 350 Pa show that for 
constant acceleration runs there is an optimal acceleration when the maximum 
cohesion force value Fc reaches its minimum (test 4). Even though the inertial force 
doubles when the acceleration increases from 0.025 m/s2 to 0.05 m/s2, the total force 
is still smaller for the higher acceleration case due to the lower cohesion force. 
Another finding is that application of linear or quadratic acceleration functions 
produce smaller maximum cohesion forces than constant acceleration. Figure 5 
shows the force from cohesion acting on the 7t boulder from regolith during the pull 
process for constant, linear, and quadratic acceleration functions. The peak value of 
the forces is lowest in the quadratic acceleration case. The force both develops and 






drops off slower than in constant acceleration case, but results in a 40% reduction in 
maximum force. 


1000 



Time (s) 

Figure 5. Force from regolith during the 8s pull tests for 7t boulder. 


Figure 5 also shows the forces with smaller cohesion strength (tests 1 and 2). We can 
expect that cohesive forces acting on a larger boulder will grow proportionally to the 
number of contacts and hence to the area of the boulder that is in contact with the 
regolith. Thus, we can estimate that the maximum force for a larger asteroid will also 
grow proportionally to its surface area and assume that 

^cl/^l =Fc2/dj 

for the same pulling method. Where Fci and Fc 2 are the maximum reaction forces 
from the regolith for the first and second boulder accordingly, and di and di are 
boulder characteristic diameters. Comparing test 4 from Table 3 to test 2 from Table 
2 for similar pulling acceleration conditions shows very good agreement with the 
above formula. Thus, the forces for larger boulders can be estimated from simulations 
of smaller boulders, which can greatly reduce the computational effort of determining 
the extraction force of large boulders. 

Conclusions 

We considered the problem of using a discrete element method simulation to estimate 
the forces that appear when pulling a boulder from an asteroid in microgravity 
conditions where cohesion forces are considered to be a significant factor. 
Simulations of 7t and 40t boulder extractions from an asteroid surface depend on the 
pulling method, cohesion strength of the regolith and boulder size. As the maximum 
force from cohesion depends on both the surface energy and geometric size of the 
grains, and we cannot model actual grain sizes due to computational limitations, we 


have to introduce an artificial surface energy depending on the particle size to match 
the desired cohesion strength of the regolith. To do this, we developed a cohesion 
strength test that uses slowly increasing artificial gravity to tear apart a regolith 
sample. The cohesion strength is then measured as a maximum stress within the 
sample during the tearing process. Using this data, we match necessary surface 
energy with cohesion strength in regolith. 

We developed a method to create a particle bed in the DEM model and perform a 
simulation of a boulder extraction in microgravity conditions. Test pulls of 7t and 40t 
boulders from an asteroid surface were performed for different boulder accelerations. 
It was found that the maximum cohesive force from the regolith depends on the 
method and speed of the pulling. For a constant acceleration pull there is the possibly 
that an optimal acceleration exists where the net maximum force from inertia and 
cohesion are smallest. Further reductions in pulling forces can be achieved by using 
linear and quadratic acceleration functions with up to a 40% reduction in force for 
quadratic acceleration. 

Comparison of extraction forces for the 7t boulder with the 40t boulder indicated that 
the extraction force for a larger boulder, using the same pulling method can be 
estimated by multiplying the smaller boulder extraction force by the ratio of square of 
characteristic diameter of the larger boulder to a square of the diameter of the smaller 
boulder: 

^cl/^l =Fc2/dj 
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